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Abstract 

The Schramm-Loewner evolution (SLE) can be simulated by dividing the time inter- 
val into N subintervals and approximating the random conformal map of the SLE by the 
composition of N random, but relatively simple, conformal maps. In the usual implemen- 
tation the time required to compute a single point on the SLE curve is O(N). We give an 
algorithm for which the time to compute a single point is 0(N P ) with p < 1. Simulations 
with k = 8/3 and k = 6 both give a value of p of approximately 0.4. 
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1 Introduction 



The Schramm-Loewner evolution (SLE) is a stochastic process that produces a random curve in 
the complex plane. In this paper we will be concerned with chordal SLE in which the random 
curve, the SLE "trace," lies in the upper half plane and goes from to oo. The classical 
Loewner equation is 

2 

9t9t(z) = —73 TT (1) 

g t (z) - U t 

with the initial condition Qq(z) = z. Here Ut is a real- valued "driving function." The equation 
defines a one parameter family g t of conformal maps from a simply connected subset of the 
upper half plane onto the upper half plane. SLE is obtained by taking U t = \fi^B t where 
k > is a parameter and B t is a Brownian motion with mean zero and variance t. For further 
discussion of SLE we refer the reader to [3, [TO] and the original references [HI Ej. 

The most common method for simulating SLE is not to numerically solve the above differ- 
ential equation. Instead one partitions the time interval into N subintervals and approximates 
gt by the composition of a sequence of N conformal maps which are approximations to the 
solution of the Loewner equation over the subintervals. Computing a point on the SLE trace 
requires evaluating the composition of roughly N conformal maps and so takes a time O(N). 
We will refer to the time it takes to compute one point on the SLE trace as the "time per 
point." 

The goal of this paper is to give an algorithm for which the time per point is 0(N P ) with 
p < 1. We do not prove that our algorithm does this. The time our algorithm takes to generate 
an SLE curve depends on the behavior of the particular curve. For certain atypical curves, our 
algorithm will not be faster than the usual algorithm. So a rigorous analysis of our algorithm 
is a daunting task. We have studied the behavior of the algorithm by simulation for k = 8/3 
and k = 6. For both of these values of k, the time per point for our algorithm is 0(N P ) with 
p around 0.4. Figures [1] and [2] show log-log plots of the time per point as a function of N for 
the usual algorithm and for our new algorithm. In the first figure k = 8/3, and in the second 
it is 6. Although the behavior of the SLE trace is qualitatively different for these two values 
of k, the behavior of the time required for the simulation is quite similar. The straight lines 
show fits by c±N for the usual algorithm and by C2N 0A for our algorithm. For iV = 100, 000 
our algorithm is faster by approximately a factor of 14 and for N = 1, 000, 000 our algorithm 
is faster by approximately a factor of 56. 

There is an inverse problem that is closely related to the present paper. Given a simple 
curve in the upper half plane, compute the corresponding driving function in the Loewner 
equation. An important application of this inverse problem is studying if a family of random 
curves is SLE. Given a large sample of the curves, one computes the corresponding samples of 
the driving function and then tests if they are a Brownian motion. This was done for interfaces 
in two-dimensional spin glass ground states in [TJ [5] and for certain isolines in two-dimensional 
models of turbulence in [21 H] . The methods of this paper apply to this inverse problem as well. 
The naive implementation of the zipper algorithm to compute the driving function for a curve 
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Figure 1: Time per point computed as a function of N, the number of subintervals in the 
partition of the time interval for k = 8/3. The top curve is the usual algorithm; the bottom 
curve is the new algorithm. The lines shown have slopes 1 and 0.4. 

with N points runs in a time 0(N 2 ). Using the methods of this paper, the algorithm runs in 
a time 0(iV 1,35 ). This fast algorithm is studied in [6]. 

In section [2J we explain the standard method for "discretizing" SLE. A particular example 
was studied by Bauer [2J. The material in this section is well known but does not seem to 
have appeared in the literature yet. In section [3] we explain our new algorithm. The algorithm 
involves several parameters, so in section H] we study how the error and time required for 
the algorithm depend on these parameters. The appendix gives some power series needed in 
particular implementations of the algorithm. C++ code implementing the new algorithm may 
be downloaded from the author's homepage. 

2 Discretizing SLE 

We are going to approximate SLE by a "discrete SLE." The term discrete SLE is a bit mis- 
leading. The random process we define produces continuous curves in the upper half plane. 
The plane is not replaced by a lattice. Let = t < t\ < t 2 < ■ ■ ■ < t n = 1 be a partition of 
the time interval [0, t]. (Thanks to the scaling property of SLE, it is no loss of generality to 
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Figure 2: Time per point computed as a function of N, the number of subintervals in the 
partition of the time interval for k = 6. The lines shown have slopes 1 and 0.4. 

take the time interval to be [0, 1].) The times tk play a special role, but the random curves are 
still defined for all time. We can think of our discrete approximation as the result of replacing 
the Brownian motion in the driving function by some stochastic process that approximates (or 
even equals) the Brownian motion at the times tfc, and is defined in between these times so 
that the Loewner equation may be solved explicitly. 

We begin by reviewing some facts about the Schramm-Loewner evolution for the upper 
half-plane, 

dt9t(z)= n 2 rn (2) 

g t {z)- yjKB t 

with go(z) = z. The set K t contains the points z in the upper half-plane for which the solution 
to this equation no longer exists at time t. Let t, s > 0. The map g t + s maps H \ K t+S onto H. 
We can do this in two stages. We first apply the map g s . This maps H \ K s onto H, and it 
maps M\K t j rS onto M\g s (K t + s \K s ). Let c/t be the conformal map that maps M\g s (K t + s \Ks) 
onto H with the usual hydrodynamic normalization. By the uniqueness of these maps, 

g s +t = gt° g s , i.e., g t = g s +t ° g^ 1 (3) 
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Then 2 2 

at dt g s+t o g s 1 (z) - U s+t g t {z) - U s+t 

Note that go(z) = z. Thus gt(z) is obtained by solving the Loewner equation with driving 
function U t = U s+t . This driving function starts at U s , and so the K t associated with g t starts 
growing at U s . In our approximation the driving function will be sufficiently nice that K t is 
just a curve which starts at U s . 

We now return to the partition — to < h < t 2 < • ■ - t n — 1, and define 

Gk = g tk o 9t k \ ( 5 ) 

So 

9t k = Gk ° G k -i o o G 2 oGi (6) 

By the remarks above, Gk(z) is obtained by solving the Loewner equation with driving function 
C/t fc _ 1+t for t — to t — Afc, where = tk —tk-i- Note that Gk maps HI minus a set "centered" 
around Ut k _ x to HI. If we consider Gk{z + U tk _ 1 +t) — U tk _ 1 , it is obtained by solving the Loewner 
equation with driving function U tk _ 1 +t — U tk _ 1 for t = to t = This driving function 

starts at and ends at 5k where 5k = Ut k — Ut k _ 1 - So this conformal map takes HI minus a set 
"centered" around the origin onto HI. 

The approximation to the SLE trace is given by ^{t) = g^iUt). Let Zk = 9t k {Ut k )- We will 
only consider the points on this curve which correspond to times t = tk- One could consider 
other points on the curve, but the distance between consecutive Zk is already of the order of 
the error in our approximation, so there is little reason to consider more points. The points Zk 
are given by 

z k = G^oG 2 1 o G- k \oG- k \U tk ) (7) 

We would like to rewrite this using conformal maps that depend only on the change in U t over 
the time intervals. Define 

h k (z)=G k \z + U tk )-U tk _ 1 (8) 

So 

z k = h 1 oh 2 o o h k -i o h k (0) (9) 

Recall that if we solve the Loewner equation with driving function U tk _ L +t — Ut k _ x for t — to 
t = Afc, we get gk(z) where 

g k (z) = G k (z + U tk _ 1 )-U tk _ 1 (10) 

Letting fk(z) = g k 1 (z), we have 

f k {z) = G k \z + U tk _ 1 )-U tk _ 1 (11) 

and so hk(z) = fk{z + 5k)- As noted above, maps HI to HI minus a curve that starts at 0. The 
driving function ends at 5k, so fk{5k) is the tip of the curve. It follows that hk(z) maps HI onto 
HI minus a curve starting at and hk(0) is the tip of the curve. Thus we have the following 
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simple picture for eq. (Q. The first map hk welds together a small interval on the real axis 
containing the origin to produce a small cut. The origin is mapped to the tip of this cut. The 
second map welds together a (possibly different) small interval in such a way that it produces 
another small cut. The original cut is moved away from the origin with its base being at the 
tip of the new cut. This process continues. Each map introduces a new small cut whose tip is 
attached to the image of the base of the previous cut. 

The key idea is to define U t for £*._i < t < tk so that gk{z) may be explicitly computed. 
There are two constraints on g^. The curve must have capacity 2 A& and must map the 
tip of the curve to 5k- Given any simple curve satisfying these two constraints and starting 
at the origin, it will be the solution of the Loewner equation for some driving function which 
goes from to 5k over the time interval [0, A&]. Different choices of this interpolating curve 
give us different discretizations. Before we discuss particular discretizations, we will discuss the 
definitions of A^ and 5k- 

The simplest choice for A^ is to use a uniform partition of the time interval, A& = 1/N. 
However, if we use a uniform partition and look at the resulting points Zk, they appear to be 
farther apart on average at the beginning of the curve than at the end. To understand this, 
consider the case of k = 8/3 which is believed to correspond to the self-avoiding walk. Let w(s) 
denote the scaling limit of the self-avoiding walk. With its natural parameterization we have 
Eus(s) 2 = cs 2u where v = 3/4 and c is a constant. If we take points at equally spaced times 
using this parameterization we will get points on the walk that are roughly equally distant. 
For SLE with its parameterization using capacity, we have Euj{t) 2 = ct for some constant c. 
To match these two parameterizations in an average sense we should take t = s 2u . Thus to get 
points approximately equally spaced on the SLE curve, we should take tk = {k/N) 2v . 

The 5k should be chosen so that the stochastic process U t will converge to y/n times Brownian 
motion as N — > oo. One possibility is take the 5k to be independent normal random variables 
with mean zero and variance nAk- If we do this, then Ut and \fnBt will have the same 
distributions if we only consider times chosen from the tk- Another possibility is to take the 
5k to be independent random variables with 5k = ii/ztA/T where the choices of + and — are 
equally probable. If we use this choice with the uniform partition of the time interval, then we 
are approximating the Brownian motion by a simple random walk. 

We now consider specific choices of the interpolating curve used for the discretization. A 
popular choice is the following. Let C be a line segment starting at the origin with a polar 
angle of air. gk maps H \ C onto H. There are two degrees of freedom for the line segment - 
its length and a. There are two constraints - the line segment must have capacity 2A& and the 
tip of the segment must get mapped to 5k- 

Consider the map 



where x, y > 0. It maps the half plane onto the half plane minus a line segment which starts at 
the origin and forms an angle a with the positive real axis. The interval [—y, x] gets mapped 
onto the slit. The length of this interval determines the length of the slit. Shifting this interval 
(relative to 0) does not change the length of the slit. To obtain the map g^, we must choose x 




(12) 
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and y so that gk satisfies the hydrodynamic normalization and has capacity 2A^. Tedious but 
straightforward calculation shows if we let 



l-Q 



'•w = (* + Wi^) ('"W^r) (13) 

then fr l {z) satisfies the hydrodynamic normalization and has capacity 2t. In particular, fk{z) is 
given by the above equation with t = We know from the general theory that gt(z) = ft(z) 
satisfies the Loewner equation ([1]) for some driving function U t - Some calculation then shows 
that 

U t = c a Vt (14) 

where 

1 - 2a 

c Q = 2^=== (15) 
ya(l — a) 

For the map gk we need U tk — U tk _ x = 8k, and so 

4 = c a \[~K k (16) 

This equation determines a. [a depends on k.) 
Define 

5 2 

v = ^- (17) 



A 



Then squaring (1161) gives 
which leads to 
and so 



16a z + va - 16a - va + 4 = (19) 



1 1 / v 

a = -±-J 20 

2 2 V 16 + v K J 



We take the choice with a < 1/2 if 8k > and the choice with a>l/2if5fc<0. Using 
hk(z) = fk(z + 8k) we find that 



l-a 



Note how confusingly similar this formula is to (fi~3l) . 
Another discretization is to take 



h k {z) = Vz 2 - 4A fc + 8 k (22) 
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This conformal map produces a vertical slit based at 5 k with capacity 2A k . It does not map the 
origin to the tip of the slit. So composing these maps does not produce a curve. Nonetheless, 
this discretization converges to SLE [2J. A vertical slit corresponds to a constant driving 
function. So this discretization corresponds to replacing the Brownian motion by a stochastic 
process that is constant on each time sub interval and jumps discontinuously at the times t k - 

3 A faster algorithm 

To motivate what we do in this section, we first consider the speed of the algorithm described 
in the previous section. Recall that points on the approximation to the SLE trace are given by 

z k = hioh 2 o o hk-i o h k (0) (23) 

The number of operations needed to compute a single z k is proportional to k. So to compute 
all the points z k with k = 1, 2, • • -N requires a time 0(N 2 ). 

It is important to note that the computation of z k does not depend on any of the other 
Zj. So we can compute a subset of the points z k if we desire. (As an extreme example, if we 
are only interested in zjy — 7(1); the time required for the computation is O(N) not 0(N 2 ).) 
For the timing tests in this paper we compute the points z^ with j — 1, 2, • • • , N/d where d is 
some integer. But we emphasize that our algorithm works for any choice of the set of points to 
compute. For the above algorithm the time grows as N 2 /d. The time per point grows as N. 
We use the time per point throughout this paper to study the efficiency. It is a natural measure 
since it depends on how finely we discretize the time interval but not on the number of points 
we choose to compute. The total time to compute the SLE trace is given by the number of 
points we want to compute on it times the time per point. Our goal is to develop an algorithm 
for which the time per point is 0(N P ) with p < 1. 

Our algorithm begins by grouping the functions in (1231) into blocks. The number of functions 
in a block will be denoted by b. Let 

H 3 = ^(i-i)6+i ° %-i)b+2 ° • • • ° h jb (24) 

If we write k as k = mb + / with < I < b, then we have 

z k = Hi o H 2 o • • • o H m o h mb+1 o h mb+2 o • • • o h mb+ i(0) (25) 

The number of compositions in (T25]) is smaller than the number in fT23|) by roughly a factor of b. 
Unfortunately, even though the hi are relatively simple, the Hj cannot be explicitly computed. 
Our strategy is to approximate the hi by functions whose compositions can be explicitly com- 
puted to give an explicit approximation to Hj. For large z, hi{z) is well approximated by its 
Laurent series about oo. One could approximate hi by truncating this Laurent series. This is 
the spirit of our approach, but our approximation is slightly different. 

Let f(z) be a conformal map from EI onto H \ 7[0, t], where 7 : [0, t] — > HI is a curve in the 
upper half plane with 7(0) = 0. We assume that /(oo) = 00, /'(oo) = 00 and /(0) = j(t). Let 
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a, b > be such that [—a, b] is mapped onto the slit 7[0,t]. So / is real valued on (—00, —a] 
and [b, 00). By the Schwartz reflection principle, / has an analytic continuation to C \ [—a, b], 
which we will simply denote by /. Let R = max{a, b}, so / is analytic on {z : \z\ > R} and 
maps 00 to itself. Thus f(l/z) is analytic on {z : < \z\ < 1/R} and our assumptions on / 
imply it has a simple pole at the origin with residue 1. So we have 



f(l/z) = l/z + J2 



c k z 



(26) 



fc=0 



This gives the Laurent series of / about 00. 



+E 

k=0 



C k Z 



(27) 



If we truncate this Laurent series, it will be a good approximation to / for large z. At first 
sight, this Laurent series is the natural approximation to use for /. However, we will use a 
different but closely related representation. 

Define f(z) = l/f(l/z). Since f(z) does not vanish on {\z\ > R}, f(z) is analytic in 
{z : \z\ < 1/R}. Our assumptions on / imply that /(0) = and /'(0) = 1. So / has a power 
series of the form 



/(*) 



2 a i 

j=0 



,z J 



(28) 



with ciq = and a± = 1. It is not hard to show that 1/R is the radius of convergence of this 
power series. We will refer to this power series as the "hat power series" of /. Note that the 
coefficients of the hat power series of / are the coefficients of the Laurent series of 1/ f ■ 

The primary advantage of this hat power series over the Laurent series is its behavior with 
respect to composition. 

{fogy{z) = l/f{l/g{z)) = f{g{z)) (29) 



Thus 



U °gY = f °g 



(30) 



Our approximation for f(z) is to replace f(z) by the truncation of its power series at order n. 
So 

/(*) ' 



/(!/*) 



.3=0 



(31) 



For each hi we compute the power series of hi to order n. We then use them and ( 1301) to 
compute the power series of Hj to order n. Let 1/Rj be the radius of convergence for the power 
series of Hj. (Rj is easy to compute. It is the smallest positive number such that Hj(Rj) and 
Hj(-Rj) are both real.) Now consider equation (|25|) . If z is large compared to Rj, then Hj(z) 
is well approximated using its hat power series. We introduce a parameter L > 1 and use the 
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hat power series to compute Hj(z) whenever \z\ > LRj. When \z\ < LRj, we just use (|24|) to 
compute Hj(z). The argument of Hj is the result of applying the previous conformal maps to 
0, and so is random. Thus whether or not we can approximate a particular Hj using its hat 
power series depends on the randomness and on which zj~ we are computing. 

4 Choosing the parameters 

Our algorithm depends on three parameters. The integer b is the number of functions in a 
block. The integer n is the order at which we truncate the hat power series. The real number 
L > 1 determines when we use the hat power series approximation to the block function. These 
three parameters control how good our approximation is. We can compute the error that arises 
from using the hat power series by computing the discretized SLE curve both using the hat 
power series and not using them. We then define the error to be the average distance between 
the points on the two curves. Given some desired level of error, we want to minimize the time 
subject to the constraint that the error is within the desired tolerance. There is a completely 
different kind of error - that introduced by approximating the Brownian motion by some other 
stochastic process. It should converge to zero as iV — > oo. The nature of this convergence and 
in particular its dependence on the method of discretization is an interesting question, but we 
do not study it in this paper. 

The behavior of our algorithm is random in that it depends on the behavior of the particular 
SLE sample. Since the behavior of SLE depends qualitatively on k, one might expect that the 
behavior of our algorithm will depend significantly on k. We have studied the algorithm for 
k = 8/3 and k = 6 and have found that the behaviors for these two values of k are remarkably 
similar. We will restrict our discussion and our plots to the case of k = 8/3, and discuss how 
K = 6 compares at the end of this section. 

We continue to use the time per point (total time divided by the number of points computed) 
as our measure of the speed of the algorithm. For this new algorithm it is essentially independent 
of d, the number of time intervals between consecutive points computed, provided d is not huge. 
The computation of the hat power series for the conformal maps Hj does not depend on how 
many points we compute. When d is large enough, the time required for this computation will 
dominate, and the time per point will no longer be independent of d. 

We first consider the effect of n. Obviously, as n grows the error should decrease and 
the time should increase. Figure [3] shows that as a function of n, the error is approximately 
proportional to L~ n . Figure H] shows the growth of the time with respect to n. (In both figures 
we have taken N = 100, 000, b = 40, and L = 4.) 

Next we consider the effect of L. As L increases we use the hat power series approximation 
only for larger z and so the error should decrease. However, using the approximation less 
frequently will increase the time required. The error decreases with L roughly as L~ n . (Figure 
[5]) The dependence of the time on L shown in figure [6] is somewhat surprising. The time does 
not grow with L as quickly as one might expect. Eventually, as L gets large the time will be 
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Figure 3: Error as a function of n, the order of the hat power series. The line is a fit by cL~ 



on the order of the time for the standard algorithm that does not use hat power series, but this 
only happens at values of L considerable larger than those shown. (In both figures we have 
taken N = 100, 000, b = 40, and n = 12.) 

Finally we consider the effect of the block size b. It is not clear a priori how the error and time 
will depend on b. Increasing b reduces the number of compositions to compute in fl25|) but it will 
also increase the radii of convergence Rj which will result in the hat power series approximation 
being used less frequently. Figure [7] shows that the error does not depend strongly on b and so 
b should just be chosen to minimize the time. This choice depends significantly on N. (In this 
figure we have taken N = 100, 000, n = 12, and L = 4.) 

Figure [8] shows the time as a function of b for three values of N. (For all three curves n = 12 
and L = 4.) The curves have similar shapes, suggesting some sort of scaling. In figure Owe plot 
the same data but now divide the time by the minimum time for that value of N and divide b 
by \/~N. The resulting three curves collapse nicely. This indicates that the optimal value of b 
is roughly proportional to \J~N . We have found that the optimal value is well approximated by 
b = 0.12VN. 

We have found that the error is typically well below L~ n . To study which value of n is 
optimal we do the following. We fix a value of n and then choose L so that L~ n = 10~ 6 . (The 
choice of 10~ 6 is ad hoc.) We then study the time per point as a function of N. Figure M 
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Figure 4: Time per point computed as a function of n 

shows the resulting plots for n = 8, 10, 12, 14. For the large values of N there is little difference 
between n = 10, 12, 14, but they are significantly better than n = 8. 

We have carried out the same simulations and generated the same plots for k — 6. Qual- 
itatively the curves are the same. The difference between the two values of k varies with the 
choices of the three parameters, but to a very crude approximation we have found that for k = 6 
the algorithm is about 20% slower, and the error is about twice as large. (The error for k = 6 
is still usually less than L~ n .) The optimal value of b for k = 6 is a bit smaller. It is better 
approximated by b = O.lv^iV. The analog of figure [TUI for k = 6 is virtually indistinguishable 
from the figure shown in which k = 8/3. In particular, the time per point is approximately 
O(N 0A ). Taking n to be 10, 12 or 14 give similar results, all significantly better than n = 8. 



A Particular hat power series 

In this appendix we give the hat power series needed for the two particular discretizations we 
have discussed. For the approximation that uses slits at an angle arc, we need to compute the 
power series of h where h(z) — (z + xi) 1 ~ a (z — x r ) a . We have 

h(z) =z{l + xizY {l - a) (1 - x r z)- a (32) 
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Figure 5: Error as a function of L, the parameter that determines how often we use the hat 
power series approximation. The line is a fit by cL~ n . 



The power series of the last two factors are given by the formula 

(i - cz)-° = £; <»(<»+ *) ••;(<»+ c k z k 



k=0 



k\ 



(33) 



For the approximation that uses vertical slits, we need to compute the hat power series of 



h(z) = \J z 1 — At + x. First consider g(z) = \J z 2 — At. We have 

1 ■ 3 • 5 • ■ ■ (2k — 1) 2k 



y/1 - Atz 2 



E 

k=0 



k\ 



(34) 



where the power series may be obtained from (|33|) with a = 1/2. Noting that h — f o g with 
f(z) = z + c, the hat power series of h is just the composition of the hat power series for / and 
g. The series for / is just 



/(*) 



1/z + c 1 + cz 



E(-D 



(35) 



m=0 
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Figure 6: Time per point computed as a function of L. 
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